#ifndef __SUBMODULAR_GRAPHCUT_H__
#define __SUBMODULAR_GRAPHCUT_H__
/*********************************************************************
 * Software License Agreement (Modified BSD License)
 *
 * Copyright (c) 2010, Daniel Munoz
 * All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 *
 *   * Redistributions of source code must retain the above copyright
 *     notice, this list of conditions and the following disclaimer.
 *   * Redistributions in binary form must reproduce the above copyright
 *     notice, this list of conditions and the following disclaimer in the
 *     documentation and/or other materials provided with the distribution.
 *   * Neither the name of the copyright holders' organizations nor the
 *     names of its contributors may be used to endorse or promote products
 *     derived from this software without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
 * HOLDERS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************/

#include <iostream>

#include <list>

// --------------------------------------------------------------
namespace submodular_graphcut
{
// --------------------------------------------------------------
/*!
 * \brief Uses graphcuts to minimize a submodular function with binary variables
 *
 * This implementation is based on a class originally written by Vladimir Kolmogorov:
 * http://www.cs.ucl.ac.uk/staff/V.Kolmogorov/software/energy-v1.1.src
 */
// --------------------------------------------------------------
class SubmodularGraphcut
{
  public:
    // for VK3, BGL
    typedef unsigned int EnergyVar;

    // for VK2
    //typedef void* EnergyVar;

    SubmodularGraphcut() :
      m_const_offset(0.0)
    {
    }

    // --------------------------------------------------------------
    /*!
     * \brief Adds a new variable to this function
     */
    // --------------------------------------------------------------
    virtual EnergyVar addVariable() = 0;

    // --------------------------------------------------------------
    /*!
     * \brief Minimizes the current energy function
     *
     * \return The minimum energy value
     */
    // --------------------------------------------------------------
    virtual double minimize() = 0;

    // --------------------------------------------------------------
    /*!
     * \brief Returns the binary value for the variable
     *
     * This should be called after calling minimize()
     *
     * \return 0 or 1 for the value of x in the minimizing solution
     */
    // --------------------------------------------------------------
    virtual unsigned int getValue(const EnergyVar& x) = 0;

    // --------------------------------------------------------------
    /*!
     * \brief Adds a unary term E(x) to this function.
     *
     * This function can be called multiple times per variable
     *
     * \param x The variable
     * \param E0 The unary energy when x = 0
     * \param E1 The unary energy when x = 1
     *
     * \return 0 on success, otherwise negative value on error
     */
    // --------------------------------------------------------------
    inline void addUnary(const EnergyVar& x,
                         double E0,
                         double E1)
    {
      // 0 = source, 1 = sink
      add_tweights(x, E1, E0); // it should be backwards
    }

    // --------------------------------------------------------------
    /*!
     * \brief Adds a submodular pairwise term E(x,y) to this function.
     *
     * This function can be called multiple times per variable pairs
     *
     * Warning, the term must be submodular:
     *         E00 + E11 <= E01 + E10
     *
     * \param x The variable
     * \param y The variable
     * \param E00 The pairwise energy when x = 0, y = 0
     * \param E01 The pairwise energy when x = 0, y = 1
     * \param E10 The pairwise energy when x = 1, y = 0
     * \param E11 The pairwise energy when x = 1, y = 1
     *
     * \return 0 on success, otherwise negative value on error
     */
    // --------------------------------------------------------------
    int addPairwise(const EnergyVar& x,
                    const EnergyVar& y,
                    double E00,
                    double E01,
                    double E10,
                    double E11);

    // --------------------------------------------------------------
    /*!
     * \brief Adds a Pn Potts high-order term E(x_1,..,x_n) over n variables, as described
     *        in Kohli et al., "P3 & Beyond: Solving Energies with Higher Order Cliques", CVPR 2007
     *
     * If implementing alpha-expansion and using 0 to represent the alpha label and 1 for the current
     * label, then Ec0 = gamma_alpha and Ec1 = gamma (where gamma = gamma_beta if all nodes in the
     * clique are currently labeled beta, and gamma = Emax otherwise), as described in Equation 38.
     *
     * \param clique_vars List of energy variables that constitute the clique
     * \param Ec0 The energy if all variables in the clique take value 0
     * \param Ec1 The energy if all variables in the clique take value 1
     * \param Emax The max energy value the clique can have
     *
     * \return 0 on success, otherwise negative value on error
     */
    // --------------------------------------------------------------
    int addPnPotts(const std::list<EnergyVar>& clique_vars,
                   double Ec0,
                   double Ec1,
                   double Emax);

    // --------------------------------------------------------------
    /*!
     * \brief Adds a Robust Pn Potts high-order term E(x_1,...,x_n) over n variables, as described
     *        in Kohli et al., "Robust Higher Order Potentials for Enforcing Label Consistency", IJCV 2009
     *
     * This method implements the alpha-expansion move for the high-order potential. It uses a basic form
     * of the potential where all nodes have equal weight for all labels and the truncation parameter Q is
     * the same for all labels.
     *
     * This method should be called when a dominant label d in the clique can be found, that is
     * when D > P - Q, where D is the number of nodes in the clique that take on label d != alpha,
     * P is the number of nodes in the clique, and Q is the truncation parameter.  Call the method
     * addRobustPottsNoDominantExpand0 if no dominant label can be found.
     *
     * IMPORTANT: this function assumes 0 represents the alpha-label and 1 represents the current
     * labeling.  Ensure to call the other energy functions appropriately.
     *
     * \param node_vars List of energy variables that constitute the clique
     * \param dominant_vars List of energy variables that currently take on the clique's dominant labels
     * \param gamma_alpha The energy if all variables in the clique take on the alpha label (value 0)
     * \param gamma_dominant The energy if dominant variables keep their label (value 1)
     * \param gamma_max The max energy value the clique can have
     * \param Q The truncation parameter, it is the number of nodes that can disagree with
     *          the dominant label.  IMPORTANT: 2Q < node_vars.size()
     *
     * \return 0 on success, otherwise negative value on error
     */
    // --------------------------------------------------------------
    int addRobustPottsDominantExpand0(const std::list<EnergyVar>& node_vars,
                                      const std::list<EnergyVar>& dominant_vars,
                                      double gamma_alpha,
                                      double gamma_dominant,
                                      double gamma_max,
                                      double Q);

    // --------------------------------------------------------------
    /*!
     * \brief Adds a Robust Pn Potts high-order term E(x_1,...,x_n) over n variables, as described
     *        in Kohli et al., "Robust Higher Order Potentials for Enforcing Label Consistency", IJCV 2009
     *
     * This method implements the alpha-expansion move for the high-order potential. It uses a basic form
     * of the potential where all nodes have equal weight for all labels and the truncation parameter Q is
     * the same for all labels.
     *
     * This method should be called when NO dominant label in the clique can be found, that is
     * there is no label d != alpha such that D > P - Q, where D is the number of nodes in the clique that take
     * on label d, P is the number of nodes in the clique, and Q is the truncation parameter.
     *
     * IMPORTANT: this function assumes 0 represents the alpha-label and 1 represents the current
     * labeling.  Ensure to call the other energy functions appropriately.
     *
     * \param node_vars List of energy variables that constitute the clique
     * \param gamma_alpha The energy if all variables in the clique take on the alpha label (value 0)
     * \param gamma_max The max energy value the clique can have
     * \param Q The truncation parameter, it is the number of nodes that can disagree with
     *          the dominant label.  IMPORTANT: 2Q < node_vars.size()
     *
     * \return 0 on success, otherwise negative value on error
     */
    // --------------------------------------------------------------
    int addRobustPottsNoDominantExpand0(const std::list<EnergyVar>& node_vars,
                                        double gamma_alpha,
                                        double gamma_max,
                                        double Q);

  protected:
    // --------------------------------------------------------------
    /*!
     * \brief Creates two edges: Source->x and x->Sink with given capacities
     *
     * Capacities can be negative.
     *
     * \param x Variable
     * \param A The capacity of edge from Source -> x
     * \param B The capacity of edge from x -> Sink
     */
    // --------------------------------------------------------------
    virtual void add_tweights(const EnergyVar& x,
                              double A,
                              double B) = 0;

    // --------------------------------------------------------------
    /*!
     * \brief Creates two edges: x->y and y->x with given capacities
     *
     * Capacities can NOT be negative.
     *
     * \param x Variable
     * \param y Variable
     * \param A The capacity of edge x->y
     * \param B The capacity of edge y->x
     */
    // --------------------------------------------------------------
    virtual void add_edge(const EnergyVar& x,
                          const EnergyVar& y,
                          double A,
                          double B) = 0;

    double m_const_offset;
};
}
#endif
